%% read rheology data
tbl = readtable('20201208_run2_CC1%_20170131_summary.xls');

%% fit model
f = @(para,X) para(1) + para(2)*X(:,1) + para(3)*X(:,2) + para(4)*X(:,1).*X(:,2);

X = [log10(tbl.ShearRate), tbl.Temperature];
Y = log10(tbl.Viscosity);

modelFit = nlinfit(X,Y,f,[1 1 1 1]);
% if exist('viscModel.log','file')
%     delete('viscModel.log');
% end
% diary('viscModel.log');
fprintf('log10(viscosity (mPa s)) = p1\n');
fprintf('                         + p2*log10(shear rate (1/s))\n');
fprintf('                         + p3*(temperature (°C))\n');
fprintf('                         + p4*log10(shear rate (1/s))*(temperature (°C)):\n');
fprintf('p1 = %g\np2 = %g\np3 = %g\np4 = %g\n',modelFit);
% diary off;

%% calculate fitted values
shearRate = logspace(2.95,4,21);
temp = 5:5:40;
viscosity = nan(length(shearRate),length(temp));
for k = 1:length(temp)
    viscosity(:,k) = 10.^f(modelFit,[log10(shearRate(:)),repmat(temp(k),length(shearRate),1)]);
end

%% plot viscosity ~ shear rate
figure(2);
hp = plot(shearRate,viscosity,'-');
hold on;
dT = 2;
for k = 1:length(temp)
    mask = (temp(k) - dT <= tbl.Temperature) & (tbl.Temperature <= temp(k) + dT);
    plot(tbl.ShearRate(mask),tbl.Viscosity(mask),'.','Color',hp(k).Color);
end
set(gca,'YScale','log');
set(gca,'XScale','log');
hold off;
legend(hp,cellstr(num2str(temp', '%g °C')),'Box','off');
xlabel('Shear rate (1/s)');
ylabel('Viscosity (mPa s)');
title('Cell carrier 1% MC');

saveas(gcf,'visc_shearRate','fig');
saveas(gcf,'visc_shearRate','png');

%% plot viscosity ~ temperature
figure(3);clf;
hp = plot(temp,viscosity,'-');
ax = gca;
ax.ColorOrder = parula(length(shearRate));

set(gca,'YScale','log');
xlabel('Temperature (°C)');
ylabel('Viscosity (mPa s)');
title('Cell carrier 1% MC');

cb = colorbar('Ticks',((0:length(shearRate)-1) + 0.5)/length(shearRate), 'TickLabels',cellstr(num2str(shearRate', '%.0f')));
colormap(parula(length(shearRate)))
cb.Label.String = 'Shear rate (1/s)';

saveas(gcf,'visc_temp','fig');
saveas(gcf,'visc_temp','png');

%% plot: relative viscosity ~ temperature normalized to 23°C
% shearRate = 964; % human
% shearRate = 526; % nyc noc
shearRate = 576; % rou aeg
temp = 5:1:40;
tempRef = 23;
viscosity = 10.^f(modelFit,[repmat(log10(shearRate),length(temp),1),temp(:)]);

figure(4);clf;
hp = plot(temp,viscosity / viscosity(temp == tempRef),'-');
set(gca,'YScale','log');
datatip(hp,'DataIndex',find(temp == 10));
datatip(hp,'DataIndex',find(temp == 23));
datatip(hp,'DataIndex',find(temp == 37));
legend(sprintf('Shear rate = %g 1/s',shearRate),'Box','off');

xlabel('Temperature (°C)');
ylabel('Relative viscosity');
title('Cell carrier 1% MC');

saveas(gcf,['relVisc_temp' '_' sprintf('%d',shearRate)],'fig');
saveas(gcf,['relVisc_temp' '_' sprintf('%d',shearRate)],'png');

%% plot: relative viscosity @ 3 temps ~ shear rate
shearRate = logspace(2.95,4,21);
temp = [10 23 37];
tempRef = 23;
viscosity = nan(length(shearRate),length(temp));
for k = 1:length(temp)
    viscosity(:,k) = 10.^f(modelFit,[log10(shearRate(:)),repmat(temp(k),length(shearRate),1)]);
end

figure(5);clf;
hp = plot(shearRate,viscosity ./ viscosity(:,temp == tempRef),'-');
set(gca,'YScale','log');
set(gca,'XScale','log');
legend(hp,cellstr(num2str(temp', '%g °C')),'Box','off');
xlabel('Shear rate (1/s)');
ylabel('Relative viscosity');
title('Cell carrier 1% MC');

saveas(gcf,'relVisc_shearRate','fig');
saveas(gcf,'relVisc_shearRate','png');
